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Abstract 



The objective of this article is to introduce the tools to analyze the con- 
trast imaging problem in Nuclear Magnetic Resonance. Optimal trajecto- 
ries can be selected among extremal solutions of the Pontryagin Maximum 
Principle applied to this Mayer type optimal problem. Such trajectories 
are associated to the question of extremizing the transfer time. Hence the 
optimal problem is reduced to the analysis of the Hamiltonian dynamics 
related to singular extremals and their optimality status. This is illus- 
trated by using the examples of cerebrospinal fluid / water and grey / 
white matter of cerebrum. 

1 Introduction 

In a series of recent articles [SJ El [7J [55] , geometric optimal control com- 
bined with adapted numerical schemes such as the Hampath code [2] is used 
to analyze the optimal control of Kossakowsky-Lindblad equations [TJ [T31 [H] . 
These equations describes the evolution of a two-level dissipative quantum sys- 
tem whose dynamics is governed by a three-dimensional system 

dx r i 

— = — 1 x + u 2 z 

at 

f t =-Ty-u x z (1) 

dz 

— = 7_ - 7+ z + my - u 2 x, 

the state variable q = (x,y,z) belonging to the Bloch ball \q\ < 1 which is 
invariant for the dynamics since the dissipative parameters A = (T,7 + ,7_) 
satisfy 2r > 7 + > |7_|. The control field is u = (ui,Ua)< The underlying 
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optimal control problem consists of minimizing the transfer time with a bound 
on the modulus of the control or of minimizing the energy transfer J \u\ 2 dt 
with a fixed control duration. 

Such a system is a model for the control of a molecule in a dissipative en- 
vironment using a laser field |22[ |2"T] but also in Nuclear Magnetic Resonance 
(NMR) spectroscopy where the dynamics of a spin 1/2 particle can be described, 
up to a renormalization, by the Bloch equation which is of the form (fTJ) with the 
restriction 7_ = 7+ |llU15ll20] . This implies that in this model, the equilibrium 
point of the free motion is the north pole (0, 0, 1) of the Bloch ball. 

In NMR, we also recall that the control is a transverse radio-frequency mag- 
netic field in the (2, y)- plane, a constant magnetic field being applied in the z- 
direction. In this domain, a striking application of geometric optimal control 
was a gain of 60 % in the control duration of the saturation of a spin 1 /2 parti- 
cle [19] . The saturation problem consists in bringing the magnetization vector 
of the sample from the equilibrium point to the center of the Bloch ball [5]. 
Such a control can be achieved by a standard NMR technique, the inversion 
recovery sequence, composed of a bang arc to invert the magnetization vector 
and a singular one along the vertical z- axis to reach the target state. It can be 
shown that the geometric time-optimal solution is the concatenation of a bang, 
a horizontal singular arc, a bang and a final vertical singular arc. The gain in 
the control duration has been shown experimentally in |19| . The experiments 
were performed using the proton spins of H2O in an organic solvent at room 
temperature. This result shows that the optimized pulse sequence can really be 
implemented with modern NMR spectrometers and a reasonable match between 
theory and experiments. 

Also this result is crucial because it confirms the ubiquity of singular trajec- 
tories in the optimal control of nonlinear systems U . In the preceding example, 
contrary to the apparent simplicity of the equations, the physical situation is 
non trivial due to the two singular directions which are necessary to compute 
the optimal solution. A direct generalization of this problem is the one of the 
contrast in NMR imaging. The model is obtained by considering two uncou- 
pled spins, each of them being solution of the Bloch equations (TT]) with different 
damping coefficients Ai = (Ti,7i), A2 = (1^2,72), but controlled by the same 
magnetic field. Denoting each system by 




= Fi(qi,Ai,u) 



where qi — (xi, yi, Zi) is the magnetization vector of each spin particle, this leads 
to a system written shortly as 



where x = ((71,92)- The associated optimal control problem is the following: 
Starting from the equilibrium point of the dynamics xq = ((0, 0, 1), (0, 0, 1)), 
the goal is to reach in a given transfer time T (which can be fixed or not), the 
final state qi(T) = for the first spin while maximizing a cost C(q2{T)) (e.g. 
\q2iT) 1 2 or the projection of q2(T) on one axis). A subcase of this problem is 
to restrict the system to X\ = x% = by considering only the component u\ of 
the control field. Our aim in this paper is to present a geometric study of this 
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control problem based on the analysis of the Hamiltonian dynamics given by 
the Pontryagin Maximum Principle (PMP) [53], the optimal control problem 
being a standard Mayer problem. 

An important point in our analysis will be the introduction of singular tra- 
jectories of the system ^ = F(x, u) whose control domain N is a smooth 
submanifold of M. p defined as follows (one can assume that N = MP): 

Definition 1 A control u £ L°°([0, T]) is called singular on [0, T] if the deriva- 
tive of the extremity mapping ere x(-) denotes 
the response to u(-) initiating from xo at t — 0, is not of full rank. 

This definition is not the standard definition in the engineering litterature, in 
particular it depends upon the control domain. But it is the correct mathemati- 
cal definition in optimal control since optimality is related to openess properties 
of the extremity mapping. 

A large amount of work has been done recently in control theory to analyze 
the role of singular extremals. This can be summarized as follows: 

1. They are feedback invariant. 

2. They can be computed using the PMP as solutions of 

. _ dH . _ dH dH 

X ~~ ^ , P ~ , U 

op ox ou 

where H(x,p,u)) = (p, F(x,u)) is the Hamiltonian lift of the system. 

As such they are extremal solutions of any Mayer type problem associated to a 
system where the cost and the boundary conditions only give boundary condi- 
tions. Also recent works have shown how to compute their first conjugate time, 
that is the first time such that the extremity mapping becomes open. This time 
corresponds also to the time where the trajectories lose their local optimality. 
Theoretically, it is related to the concept of singularity of Lagrangian manifolds 
[3] and is numerically implemented in the Hampath code |14| . 

Hence going back to the contrast imaging problem, a research program is 
to analyze the Hamiltonian dynamics of the singular extremals completed by 
numerical simulations to compute the optimal solutions. This is a difficult 
task since the problem is depending upon different relaxation parameters in the 
Bloch equation. In this paper, we will present the geometric tools and some 
preliminary numerical results in two particular cases by considering only one 
component of the control field. 

The organization of this article is the following. In the first section, the 
Maximum Principle is introduced to select minimizers among extremal solutions 
in a Mayer problem. The role of singular extremals is presented and their 
optimality status is determined using the concept of conjugate points. In a 
second section, a thorough analysis of the geometric control of a single spin 
1/2 particle is presented and it plays for specific values of the parameters, an 
important role in the problem. In the final section, we numerically analyze 
the geometry of singular extremals in view of studying some specific cases in 
NMR. Numerical computations of the optimal solution are also presented for 
two regularized cost functionals. 
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2 Geometric optimal control 



2.1 Preliminaries 

One considers a Mayer problem given by the following data : 

1. A smooth system ^§ = F(x,u), x £ W 1 with fixed initial state Xq and a 
transfer time T, the controls being the set U = L°°([0, T], U) of bounded 
measurable mappings valued in a control domain U C M p . 

2. A terminal manifold M defined by f(x) = where / : E n — > is a 
smooth mapping. 

3. A cost to minimize : min u (.) e ^ C(q(T)) where C : R" — > IR is a smooth 
regular mapping. 

The geometric setting is the following. Denote x(t, u) the trajectory initiating 
from Xq and associated to u, A(xq,T) = \J u ^ux(T,u) the accessibility set at 
time T and introducing the manifold C m = {/ = 0, C(x) — m} where m is a 
parameter, an optimal control u* is such that x*(T) = x(T,u*) belongs to the 
boundary of A(xq,T), f(x*(T,u*)) = and m is minimum. 



2.2 Pontryagin Maximum Principle 

The application of the maximum principle leads to the following necessary con- 
ditions |23j . 

Proposition 1 Let u*(-) be an admissible control whose corresponding trajec- 
tory x*(t) = x(t,u*) is optimal. Then there exists an absolutely continuous 
vector function p*(-) and a scalar po < such that if we denote by H the 
pseudo-Hamiltonian H(x,p, u) = (p, F(x, u)) , the following necessary conditions 
are satisfied a.e. on [0,T]: 

dx* dH , , dp* dH . , 

H(x*,p*, u*) = m&x H (x* , p* , it) (3) 

together with the boundary conditions: 

f(x*(T)) = (4) 
P*(T) = Po ^(x*(T)) + <£, %{x*{T))), (5) 

£ £ K fe , po < (transversality conditions) . 

Definition 2 We call extremals a triplet (x,p, u) solution of HP and of the max- 
imization condition It is called a BC- extremal if it satisfies the boundary 
conditions and |3p. 
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2.3 A review of the properties of singular trajectories 

Next we present some concepts and properties about singular trajectories which 
are important in our analysis, see [1] for a complete presentation. 

We have the following characterization of singular control which allows a 
practical computation. 

Proposition 2 The control u(-) and the corresponding trajectory x(-) are sin- 
gular on [0, T] if and only if there exists a non zero adjoint vector p(-) such that 
(x,p,u) is solution a.e. on [0, T] of 

x-— p--— — -0 (6) 

dp ' dx ' du 

where H(x,p,u) — (p, F(x,u)) is the Hamiltonian lift. Moreover for each < 
t <T, p(t) is orthogonal to ImE' Xo,t (u\[Q tt ])- 

Definition 3 A singular extremal is a triple (x,p,u) solution of the above equa- 
tions. It is called: 

1 . Regular if 0^ is of maximal rank. 

2. Strongly normal if for each <t\ <t 2 <T, I m E' x( - tl ^ ta ~ tl H[ tl ,t 2 ]) is of 
corank one. 

3. Exceptional if H = 0. 

Computation in the regular case: Using the condition ^ 0, one can 
solve locally the equation ^ = and compute the singular control as a function 
u(z), z = (x,p) and plugging such u in H defines a true Hamiltonian denoted 
again H(z). If II is the standard projection (x,p) M> x, one can define the 
exponential mapping exp^ : {t,p) h4 H(exp[tH(xo,p)]) where xq is fixed. This 
leads to the following definition. 

Definition 4 Let z(t) = (x(t),p(t)) be the reference extremal solution of H . 
The time t c is said to be geometrically conjugate if exjp Xg is not of maximal 
rank at (t c ,p(0)). 

We have the following standard test: 

Proposition 3 The time t c is geometrically conjugate if and only if there ex- 
ists a non trivial Jacobi field J(t) solution of the variational equation 5z — 
dH(z(t))Sz and vertical at time and t c : dTL(J(0)) = dH(J(t c )) = 0. 

The following result is crucial in our optimality analysis: 

Proposition 4 In the strongly normal case and in the non exceptional situa- 
tion, the extremity mapping E X °' T is open for the L°° - topology at W|[o,t] where 
t > ti c . 

Application: One consider a control system of the form F(x,u) — Fq(x) + 
u\F\{x) + u%F2{x) where the control domain U is the disk u\ + u\ < 1. The 
Hamiltonian is H = H0+U1H1+U2H2 where Hi = (p, Fi(x)). The maximization 
condition © leads to 

^r*- 1 ' 2 ( 7 ) 
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outside the switching surface S: H\ = H2 = 0. The corresponding extremals 
are called of order zero and there are solutions of the smooth vector field defined 
by H{z) = H + \J H\ + iff. The corresponding solutions are regular singular 
extremals if one restricts the control domain to the unit sphere S 1 . Introducing 
u\ = cos a, U2 = sin a and extending the system using a = v, they correspond 
to singular trajectories of the extended system: 

x = Fq + cos aF\ + sin olF-i , a = v. 

The case of afflne systems: For optimality analysis, one restricts our study 
to a single input affine system: x = F + u\Fi, \m\ < 1. Relaxing the control 
bound, singular trajectories are parameterized by the constrained Hamiltonian 
system: 

. dH . OH dH TT 
x = -s-, v = --5-, ^— =Hx = 0. 
dp ox ou\ 

The singular extremals are not regular and the constraint H\ = has to be 
differentiated along an extremal to compute the controls. Introducing the Lie 
brackets of two vector fields X, Y computed with the convention 

f)X BY 
\ X ^]{x) = ^{x)Y{x)- -^(x)X(x), 

and related to the Poisson bracket of the Hamiltonian lifts Hx , Hy by the rule 
{H X ,H Y } = H [X ,Y], one gets: 

H, = {Hi, H ] - {{H 1 ,H },Ho} + tii{{H 1 ,H f) },H 1 } = 0. 

A singular extremal such that {{Hi, H }, Hi} ^ is called of minimal order 
and the corresponding control is given by 

{{Hi, H }, H } 
Uls -~{{Hi,H },Hiy (8) 

Plugging such ui s into H defined a true Hamiltonian, whose solutions initiating 
from Hi = {Hi,H } = defined the singular extremals of order zero. They are 
related to the regular case using the following Goh transformation. Assuming 
Fi non zero, then there exists a coordinate system (xi,X2, ■ • ■ ,x n ) on an open 
set V such that Fi = and the system splits into: 

x = F'(x',x n ), x n = Fq(x') +ui 

where x' — (xi, ■ ■ ■ , x n -i) and the system F' defined on an open subset V where 
x n is taken as the control variable is called the reduced system. We introduce 
the reduced Hamiltonian H'(x',p', x n ) — (p' , F'(x', x n )). One has: 

ddH -fH H\- 9H ' m 
di^- {Hl > H ° } --^x- n (9) 

= m>*oh*} = -^ do) 

This gives the relation between the affine singular case and the regular one. 
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2.4 High-order maximum principle in the affine case 



As a consequence and using the generalized Legendre Clebsch condition deduced 
from the high-order maximum principle |17j . one gets the following. 

Consider the Mayer problem for an affine system of the form x = Fq(x) + 
u\F\{x), \u\ \ < 1. Then the following conditions are necessary for optimality: 

, _ dH , _ dH 

op ox 

H(x, p, u) = max H (x, p, v) 
\v\<i 

with the boundary conditions 

f(x(T)) = 0, p(T) = p ^- + (£, g>, po < 0. 

Moreover if the control is singular and non saturating, i.e. |ui s | < 1, the gener- 
alized Legendre-Clebsch condition must hold: 

{{Hi,h },Hi}>o. (11) 



2.5 Generic classification of the bang-bang extremals near 
the switching surface 

An important issue in the contrast problem is to apply the results from [18] to 
classify the extremal curves near the switching surface. The switching surface 
is the set E : Hi = 0, while the switching function is t H> $>(z(i)) = Hi(z(t)), 
where z(t) is an extremal curve. Let E s : Hi = = {Hi, Ho}- The singular 
extremals are entirely contained in E s . A bang-bang extremal z(t) on [0, T] is an 
extremal curve with a finite number of switching times < ii < • • • < t n < T. 
We denote by £ + , £_ the regular arcs for which u = ±1 and by £ s a singular 
arc; £i£2 denotes an arc £i followed by an arc in- 
ordinary switching time. It is a time t such that a bang-bang arc switches 
with the condition $(t) = and 4>(i) = {Hi,H Q } ^ 0. According to the 
maximum principle near E, the extremal is of the form if $(t) > and 

if $(f) < 0. 

Fold case. It is the case where a bang arc has a contact of order 2 with the 
switching surface. Denoting $± = {{Hi, Ho}, Hq} ± {{Hi, Ho}, Hi} the second 
derivative of the switching function, if non zero, we have three cases: 

1. Hyperbolic case: At the switching point, one has > and $_ < 0. 
At E s , a connection is possible with a singular extremal which is strictly 
admissible and satisfies the strong Legendre-Clebsch condition. The ex- 
tremals are bang-singular-bang £±£s£±- 

2. Elliptic case: At the switching point, one has < and $_ > 0. A con- 
nection with the singular extremal is not possible and every extremal curve 
is bang-bang but with no uniform bound on the number of switchings. 

3. Parabolic case: It is the situation where and $_ have the same sign 
at the switching point. One can check that the singular extremal is not 
admissible and every extremal curve near the switching point is bang-bang 
with at most two switchings, i.e. or 
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2.6 The concept of conjugate points in the affine case 



According to [I] , this concept is related to the notion of conjugate points in the 
regular case using the Goh reduction. The important property is the following 
geometric characterization. Let z(-) — (x(-),p(-)) be a singular extremal associ- 
ated to the control defined by ©. Assuming that it is strictly admissible, one 
can embed the singular extremal into a surface S formed by all the singular 
extremals starting from xq = II(z(0)) and with initial adjoint vector p such 
that \p — p(0)\ < s. Up to the first conjugate point, the extremal synthesis is 
bang-singular-bang £±£ s £±, where bang arcs will be in the neighborhood of the 
reference singular extremal, which is related to the problem of extremizing the 
transfer time and hence to the Mayer problem. This synthesis is also valid in a 
C°- neighborhood of the reference case in the limit situation where m — > +oo, 
m being the control bound. 

2.7 Application to the contrast problem 

A direct application is the contrast problem with the boundary condition q\ (T) = 
and the cost \q2(T)\ 2 . Splitting the adjoint vector into p — {p\,P2), we deduce 
the transversality condition pz(T) = —2p q 2 {T), p n < 0. The case p — gives 
P2(T) = 0. Since the system splits into: 



where p = (pi,P2) is written as a row vector. The condition P2(T) = corre- 
sponds to a second spin which is not controlled. In the non trivial case, po is 
non zero and it can be normalized to po = —1/2. 

2.8 The embedding results 

From the previous results, one deduces the following propositions. 

Proposition 5 The time minimizing solutions of the first spin 1/2 particle can 
be embedded as extremals of the contrast problem, with p$ = 0. 

Proposition 6 In the contrast problem, the extremals of the single-input case 
are extremals of the bi-input case. 

3 The single spin 1/2 case 

Since in the contrast problem, the magnetization vector of the first particle 
has to be set to 0, an important issue is to analyze this task and the underlying 
problem of reaching this target in minimum time. Besides, the optimal solutions 
of such a problem can be embedded into the extremal solutions of the contrast 
problem. Indeed, if the transfer time in the contrast problem is exactly this 
minimum time, they are the only solutions satisfying the boundary conditions. 
Hence, in this section, based on the preliminary work [TJJ], we make a thorough 
analysis of the single input case, with an emphasis put on the role of singular 
trajectories. 



91 = F i{qi, u ), Q2 

the adjoint system decomposes into: 



^2(<Z2, u) 





3.1 Preliminaries 



First of all, since the initial condition is on the z- axis of revolution of the 
system, the control problem can be restricted to the 2D- meridian of the Bloch 
ball and the control field reduced to only one component [SJ[7]. The system is 
q = F (q) +u 1 Fi(q), \m\ < to, where 

d d 
Fo = -Tyj- + 7 (1-*) — 
oy Oz 

d d 
Fi = -z— + Vjr- 
oy oz 

Denoting S = 7 — T, the following Lie brackets are relevant in our analysis: 
[ F 1 ,F } = (- 1 + 5z)-^ + Sy^- 
[[FuFolFo] = (7(7 - 2r) - S 2 z)^-+S 2 y-^ 
[[F 1 ,F ],F 1 ]=2Sy— + ( 1 -25z) — . 

3.2 Singular trajectories and optimality 

The singular trajectories are located on the set S : det(Fi, [Fi,Fo]) = 0, which 
is given in our case by y(—2Sz + 7) =0. Hence it is formed by the z- axis 
of revolution y = and the horizontal direction z = "f/(2S). The singular 
control is given by D' + u u D = where D = det( J Fi, [[F l: Fo], Ft]) and D' = 
detiFuHF^FolFo}). 

• For y = 0, one has D ~ —2(7 — 26z) and D' = 0. The singular control is 
zero and the singular arc is solution of 

y = -y, z = t(! - z ) 

where the equilibrium point (0,1) is stable if 7 7^ 0. 

• For z = 7/(25), D = -26y 2 , D 1 = y 1 {2T- 1 ) and u ls = 7 (2T - 7)/(2<5y), 
2r — 7 > 0. Hence along the horizontal direction, the flow is 

. 7 2 (2r- 7 ) 

y = ~ T y- A5 2 y ' 

and |wi s | — > +00 when y — >• 0. 

More precisely, along the horizontal singular line, the following proposition is 
crucial. 

Proposition 7 Ifj^O, the singular control along the horizontal singular line 
is in the L 1 but not in the I?- category, near y = 0. 

This can be straightforwardly shown by using the relations: 



[ u ls (t) 2 dt= [ u ls {y)^- (12) 

Jt Jy V 
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where to and t\ are the initial and final times along the singular arc and yo the 
initial y- coordinate of this arc. One deduces that the integrand of Eq. ([12")) 
scales as 1 jy when y —> and that the corresponding integral has a logarithmic 
divergence. 

In order to study the optimality of the singular directions, one uses the 
generalized Legendre-Clebsch condition, which takes the following form for a 
2D-system. Let D" = det(Fi,F ) = jz{z - 1) + Ty 2 . The set C = {D" = 0} is 
the collinear set. If 7 ^ 0, this set is not reduced to a point and the intersection 
with the horizontal singular line is empty, except in the case 7 = 2T . Singular 
lines are fast if DD" > and slow if DD" < 0. 

To complete the optimality analysis, one introduces the clock form lj = pdq 
which is defined outside the collinearity set C by the relations (p,Fq) = 1 
and (p,Fi) = 0, the sign of dia being given by 2/(7 — 25z).This form allows to 
deduce the optimality of singular extremals and to compare two different regular 
extremals when they do not cross the singular and collinearity sets [3]. 
Parameters conditions 

The interesting case is when the horizontal singular line zq = 7/(25) cuts the 
Bloch ball |q| < 1, which gives the condition T > 3-f/2 and — 1 < zq < 0. Using 
the generalized Legendre-Clebsch condition, one deduces that the horizontal line 
is optimal and the z- axis of revolution is optimal if 1 > z > zq. In particular, 
this line is slow in the domain zq > z > — 1. Using the clock form, one can 
deduce that near the origin, the broken singular arc formed by a horizontal arc 
followed by a vertical line is time-minimal for the unbounded case, provided 
admissible controls are extended to L 1 . Note also that such a broken singular 
trajectory is not in L 2 and is not optimal for the energy minimization problem 
®. 

Having made this optimality analysis, one can deduce the time minimal 
optimal synthesis near the origin which is introduced next. 

3.3 The SiSi singularity (Interaction between two singular 
arcs) 

Assume 7 ^ 0, F > |7, < m and the control bound is large enough such 
that the bang arc u\ = m starting from the north pole (0, 1) intersects the 
horizontal singular arc zq — 7/ (25) at a point A. The horizontal singular line is 
admissible up to a saturating point B. The time minimal synthesis, with initial 
point A, is represented on Fig. [TJ Due to the saturation phenomenon at B, there 
is a birth of a switching locus E3, but the remarkable fact due to the interaction 
between the horizontal and the vertical fast singular directions is the following 
concept. 

Definition 5 We call bridge between the horizontal singular arc and the vertical 
singular one, the bang arc, such that the concatenation singular-bang- singular 
is optimal. 

This concept is important and leads to a generalization in higher dimension, 
which plays an important role in the contrast problem. 

In order to compute the global optimal synthesis provided m is large enough, 
we must analyze the synthesis near the north pole, which is presented next. 
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Figure 1: Schematic representation of the optimal synthesis when the initial 
point of the dynamics is the north pole. An arbitrary zoom has been used to 
construct the figure. Regular curves are plotted in blue (dark) and red (dark 
grey) for control fields equal to — m and +to, respectively. The optimal singular 
trajectories are displayed in green (light grey). The black line is the switching 
curve, while the dashed one is the non-admissible part of the horizontal singular 
line. 

3.4 The SiCo singularity (Interaction between the collinear 
set and the singular set) 

Observe that the north pole is a stable fixed point for the free motion and the 
vertical singular direction is a fast direction, near the north pole, provided z < 1, 
as a consequence of the strong Legendre-Clebsch condition. The collinear set C 
corresponds to an oval below the line z = 1. Using polar coordinates y = r sin cj), 
z = r cos <f>, one gets 



Hence, r which represents the purity of the quantum system decreases outside 
the oval and increases inside. The north pole is a singularity of the optimal prob- 
lem which combines a collinear situation with a singular one, but the analysis 
of the time-minimal synthesis near this point is simple because of the symmetry 
of revolution. Indeed, the only way to leave this singularity is to use a bang arc 
u = ±m, which gives a boundary arc of the accessibility set. This first bang arc 
is followed by another bang arc u — =pm to fill the interior of the domain and 
to reach the vertical singular axis. This gives an optimal bang-bang policy (see 
the top part of Fig. [IJ. 

3.5 The global synthesis 

Under our assumptions (r > 37/2 > 0, to large enough), the global time min- 
imal synthesis starting from the north pole is easily obtained gluing the two 



dr 



Ty 2 - 7(2 



1) = -D". 
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previous syntheses (the one associated to the SiSi case with the one of the SiCo 
case). It is represented on Fig. [T] The switching locus is formed by the arc 
starting from the north pole and reaching the horizontal singular arc at A (it 
is denoted Si in the figure), the horizontal singular segment £2 between the 
points A and B, the switching locus £3 due to the saturation phenomenon and 
the part of the vertical singular direction between D and (the £4 segment), 
D being the extremity of the bridge. The bang arc with u = —m starting from 
A is separating the two domains, one with a bang-bang policy and the other 
containing a non trivial singular arc. 

At the limit, when m — > +00, it gives the synthesis constructed in Ref. [19 
where the total time to reach the origin is formed by the time to follow the 
broken singular-singular arc between A and 0. Observe also that according to 
our analysis, the usual policy in NMR, the inversion recovery sequence, where 
only the vertical singular arc is used, is slow if z < zq. 

Also, note that the switching locus has a complicated structure, but due to 
the symmetry of revolution, all the cut points, i.e. the first points where the 
extremal trajectories cease to be optimal, are on the vertical z- axis where two 
symmetric solutions starting respectively on the left and right part of the Bloch 
disk intersect at the same time. 



4 Preliminary results in the contrast problem 

As mentioned in the introduction, the goal of the contrast problem is to bring 
the magnetization vector of spin 1 towards the center of the Bloch ball together 
maximizing the modulus of the magnetization vector of the other specie. Note 
that such a computation could have potential applications in magnetic resonance 
imaging in order to optimize the contrast of a given imaging [8, 9J. Roughly 
speaking, the species with a zero magnetization will appear dark, while the other 
species with a maximum modulus of the magnetization vector will be white. We 
introduce in the following a simple model reproducing the main features of this 
control problem. We describe the general structure of the optimal solution and 
we compute them for two particular examples. 



4.1 The model system 

Each spin 1/2 particle is governed by the Bloch equation: 

dM x M-r 



dt T 5 ' v 



+ u v M z 

2 



dt T 2 
dM z (M -M z ) 

—7— = = 1" U x My - UJyM x 

at J 1 

where the state variable is the magnetization vector and T\ , Ti are the relaxation 
times. The control is the magnetic field uj = (lo x , uj y ) which is bounded by < 
Umax- We use the normalization introduced in [19J. The normalized coordinates 
are q = (x,y,z) — (M x ,M y ,M z )/Mo. In these coordinates, the equilibrium 
point is the north pole (0,0,1) and the normalized control is u — (u x ,u y ) — 
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^ 27r (bj x ,ujy), \u\ < 2n, while the normalized time is given by r = ui max t / (2tt) . 
Hence the system takes the form: 

x = —Tx + UyZ 

y = -Ty - u x z 

z = 7(1 - z) + (u x y - u y x) 

where T = 2ix j {io max T2) and 7 = 2*k / (u) max Ti). In the experiments, Lo max can 
be chosen up to 15 000 Hz but the value 2tt x 32.3 Hz will be considered in 
this paper. The experiments are done for the contrast problems of the cere- 
brospinal fluid/ water [25] and the grey/ white matter of cerebrum cases [TU]. In 
the cerebrospinal fluid/water situation, the relaxation parameters for the first 
spin describing the fluid are T\ = 2000 ms and T2 = 200 ms, while for the 
second spin T\ = T2 = 2500 ms. In the second example, the rates of the grey 
matter are taken to be T\ = 920 ms and T2 = 100 ms, the rates for the white 
matter being T\ = 780 ms and T2 = 90 ms. 

4.2 Computation of the singular flow 

One restricts to the situation where the control field has only one component and 
the contrast problem is governed by the differential system q — Fo(q) + uiFi(q), 

q = (yi,zi,y2,z 2 )- 

E 2 , d d . 

Denoting Si = Ji — i = 1, 2 one has: 

[F 1 ,F } = ^(- 7i + S i z i )^ + Siy i £- 

[[F U F , Fo}] = X)[7*(7i - 2r<) - + Sfy~ 

[[F 1 ,F },F ] = 2S ^dy- + (7i - 2S ^d^ 

and the corresponding singular flow is defined by: 

E x = {H!,H } = {{H U H }, H ] + u ls {{H u H }, H x ) = 0. 

Since the equations are linear with respect to p, for each initial condition 
qo, this defines a two-dimensional surface S(qo) in the state space. An ad- 
ditional condition is provided by the generalized Legendre-Clebsch condition: 
{{Hi, Ho}, Hi} > 0. The structure of this surface is related to the relaxation 
parameters (rj,7j). 

If the transfer time is not fixed, this leads to the additional constraints 
Hq = 0. In this case, the singular flow defines a single vector field in the state 
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space, since the adjoint vector can be eliminated and the restricted singular 
control is given by: 

D'(q) 

where 

D' (q) = det(F , F 1 , [Ft , F ] ,[[F lt F ] , F ] ) 
D(q) = det(F , F t , [Fi,F„], [[F U F ], Fi]) 

with the corresponding vector field 

^ = Fo(<z) -^y Fl(9) 

which can be analyzed using the time reparameterization dr = dt/ D{q{r)). In 
this framework, singular trajectories are used to classify the systems. 

In the general similar computation shows that the singular trajecto- 

ries are solutions of an equation of the form: 

M = Fo{q) - ^(a-X) Fliq) (13) 

where A is a one dimensional time dependent parameter whose dynamics is 
deduced from the adjoint equation. The solutions of Eq. (|13[) emanating from 
<7o will form S(qo). 



4.3 Numerical simulations on singular trajectories 

We present some numerical simulations concerning the singular trajectories. 
The projection of S(qo) on the planes (2/1,21), (2/2,^2) shows the effect of the 
relaxation parameters on the contrast. This point is illustrated by the figures [5] 
and[3]for the cerebrospinal fluid/water and grey/white matter of cerebrum cases, 
respectively. In each example, we assume that a bang pulse of large amplitude 
has been first applied to the system, the initial point of the singular flow is 
of coordinates ((— y/\ — Zq, zq), (— yl — Zq, Zo)) where z — zq is the horizontal 
singular line of the first spin. This first bang is necessary so that the singular 
trajectory of the spin 1 can reach the center of the Bloch ball. One clearly sees 
in Fig. [3] the similar structure of the different singular trajectories of the two 
spins. The situation is completely different in Fig. [5] for the first example. This 
explains the excellent and weak contrasts that can be reached in the first and 
second examples with an optimal sequence of the form bang-singular. Note that 
some singular control fields diverge as displayed in Figs. [5]and[21 The conjugate 
points defined in Sec. [5] have been computed for each singular extremal as 
shown in Fig. 3) Similar results have been obtained for the spin 2 and for 
the grey /white matter case. This shows that the structure bang-singular is 
not optimal since the first conjugate point occurs before the saturation of the 
spin. A more complicated pulse sequence such as bang-singular-bang-singular 
has therefore to be used. 
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Figure 2: Structure of the projection of the singular flow onto the planes (yi, zi) 
and (7/2, Z2) in the cerebrospinal fluid/water case. The trajectories are plotted 
in black (solid line) and in red (dashed line). The control fields of the dashed 
extremals diverge. The trajectories have been plotted up to the explosion of the 
field (The absolute value of the field is larger than 10 5 ). The horizontal solid 
line is a singular line of the first spin. 



15 



Figure 3: Same as Fig. [2] but for the grey/white matter case. An explosion of 
the control field is observed for the red trajectories. 



HI 
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Figure 4: Zoom of the results of Fig. [2] for the spin 1 near the origin. The 
red crosses indicate the position of the first conjugate point. The dashed lines 
represent the singular trajectories for which the control field diverges. 
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4.4 Some preliminaries numerical results on the contrast 
problem 

Due to the numerical difficulty of the computation of the bang-singular-bang- 
singular optimal sequence, we only present in this section some preliminary 
results. To make the numerical simulations, we have used a differential con- 
tinuation method of the Hampath code [14] where the cost is regularized by 
adding a L 2 (or a L 2 ~ A ) penalty on the control. Note that another continuation 
on the transfer time has also been used in the computations. Such results can 
be compared with the GRAPE algorithm [T^l [TH1 HH [25] which is a standard 
approach in NMR to solve the optimization problems. 
In the first study, the cost is regularized as 



where A is a continuation parameter and the transfer time varies starting from 
T m in + e to 2T m j„, where T m i n is the minimum time to saturate the first spin 
and £ < 1 an arbitrary parameter. According to Sec. [3J in the limit case where 
T = T m i n > the optimal solution of the contrast problem is exactly the solution 
of driving the first spin to the origin. 

The different numerical results are presented in Figs. [3] El and [5] The dif- 
ferent behaviors for the two examples can be clearly seen since the best contrast 
V 2/1(0 + z ! W 1S °f the order of 0.73 and 0.07 in the first and second examples, 
respectively. Note also that for T = T m . m + e, the trajectory of the spin 1 is 
very close to the trajectory for saturating this spin in minimum time. 

An interesting phenomenon can be observed in the cerebrospinal fluid/water 
situation in Fig. where there exists a bifurcation of the optimal policy when 
T increases. This is related to the introduction of a bang-bang policy associated 
to a SiCo singularity. Also we observe that the optimal policy is crossing the 
Z\- axis of revolution. Further work is necessary to improve the continuation 
method near A = 1 since the L 2 - regularization of the cost is not adapted to the 
control saturation that can be found in the SiSi singularity where the control is 
L 1 and not L 2 . 

This point is illustrated by a second series of simulations where we have 
considered the following regularized cost: 



and as before a continuation has been performed on the control duration T. The 
computation has been done for A = 0.9. Similar contrasts have been reached 
in this second situation. Note, however, the different peaks appearing in the 
evolution of u, to be compared to the first regularization. 

We complete this paper by illustrating our numerical results on a simulated 
and a real contrast experiments. For the simulated experiment, we consider two 
surfaces as displayed in Fig. [TJJ filled in with spins 1 or 2 in a homogeneous 
manner. We apply the optimal control field and we associate a color to the final 
modulus of the magnetization vector of the spin 2. This color is white if the 
modulus is equal to 1, black if it is zero and a grey variant between. One clearly 
sees in Fig. [TTJ the excellent and weak contrasts that can be obtained in the 
first and second examples. 
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Figure 5: (The cerebrospinal fluid/water case) Trajectories of the first and sec- 
ond spins for T m i„ + e and 2T TO ,- n in solid and dashed lines, respectively. The 
horizontal singular line is plotted in solid line. The parameter A is taken as 0.9. 
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Figure 6: (top) Evolution of the control field associated to Fig. \5\for T m i n +e and 
2T m i n in solid and dashed lines, respectively. The time T has been normalized 
to 1 to plot the two control fields on the same figure, (bottom) Evolution of the 
contrast parameter y/ i/2(T) 2 + Z2(T) 2 as a function of the control duration. 
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Figure 7: Same as Fig. [3]but for the grey/white matter of cerebrum. 
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Figure 9: Same as Fig. [5] but for the second regularized cost functional. The 
parameter A is taken as 0.93. 
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Figure 11: Simulated experimental results on the contrast problems of the cere- 
brospinal fluid/water (middle) and the grey/white matter of cerebrum (bottom) 
examples. The inner disk mimics the spin 1, while the outside ring mimics the 
spin 2. The two surfaces are separated by a thin black circle. The top figure 
is a reference image where a 90 degree pulse has been applied to the two spins. 
The middle and bottom images are a representation of the contrast as could be 
done in a real experiment. For these two images, the control sequence is the 
optimal field. A color has been associated to each value of the contrast between 
and 1, and 1 corresponding respectively to the colors black and white. 
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In Fig. 1121 we compute by interpolation the contrast between the two spin 
particles for different values of the relaxation parameters. To have a 2-D rep- 
resentation, one fixes the first spin parameters. In the top figure, the first spin 
corresponds to the cerebrospinal fluid (If — 2000 ms and T\ — 200ms), in 
the middle figure, it is the gray matter (T\ — 920 ms and T\ = 100 ms) and 
in the bottom one, it is the deoxygenated blood (T-j 1 = 1350 ms and T\ = 50 
ms), this latter example being illustrated below experimentally. In each case, 
we fix the control duration T to 1.5T m i n and we choose the regularized cost 
C{x(T)) + (1 - A) J Q T \u\ 2 - x {t)dt, with A = 0.9. We consider the following vari- 
ations for the parameters of the second spin: 

Umin ^2 — Vmaxi 

where (x min , x max ,y mm , y max ) = (80, 4000, 160, 4000) for the fluid case, (45, 1500, 90, 1500) 
for the matter case and (20,2000,40,2000) for the blood case. The linear in- 
equalities T2 < 2Ti due to the physical model, leads to convex polyhedrons 
in Fig. [T2j The starting point of the forthcoming homotopies corresponds to 
S = (Tf, T^) for which the contrast is zero. Then we discretize the edges of the 
polytope P into n points and for each = (if' , T^). k = 1, n, we perform 
a linear homotopy from S to i*& by introducing a parameter A such that : 

if = Tl + \(T?' k - Tl) 
T% = Ti + \(T^ k -Ti). 

At the end, we have n lines starting from S which mesh P and to complete the 
figures, we use a standard Matlab interpolation function. We can see in Fig. 
IT21 on the top figure, that the contrast in O = (2500, 2500), which corresponds 
to the Fluid/Water case, is nearly 0.7 and on the middle one, in O = (780, 90), 
the Gray/White matter case, the contrast is almost 0.1, which agree with the 
results given in Figs. ITD1 and [51 respectively. 

4.5 Some preliminary experimental results 

The first experimental results on the contrast problem are represented on Fig. Q2] 
and correspond to samples reproducing the case of the deoxygenated/oxygenated 
blood. Such results can be compared to Fig. [TT] where they have been numer- 
ically simulated for other samples in an ideal experiment. The preliminary 
experimental results are promising even if some artefacts due to the inhomo- 
geneities of the magnetic field deteriorate the quality of the image. 

5 Conclusion 

In the conclusion, we discuss some important issues related to our study. 

Mathematical problems. The important remaining question is to analyze 
the dynamics of the singular flow in relation with the relaxation times and in 
particular the asymptotic of the trajectories. 
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Figure 12: Interpolation results for the (top) Fluid case, the (middle) Gray 
matter case and the (bottom) Blood case. The contrast is computed with respect 
to the spin 2 relaxation parameters. The parameters of the spin 1 are fixed to 
the ones of the points S. The points O give us the contrast for known problems, 
as the fluid/water case on the top figure, gray/white matter case on the middle 
one and deoxygenated/oxygenated blood on the bottom figure. 
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Figure 13: Experimental results: The inner circle shape sample mimics the 
deoxygenated blood, where T\ = 1.3 s and T2 = 50 ms; the outside moon shape 
sample corresponds to the oxygenated blood, where T\ = 1.3 s and T 2 = 200 
ms. The goal of the control is to saturate the inner sample and to maximize the 
remaining magnetization of the outside sample. The upper image is a reference 
image after a short 90 degree pulse on both samples. The image at the middle is 
the remaining Y magnetization \M y \ after the optimized pulse, the lower image 
is the remaining Z magnetization \M Z after the optimized pulse. 
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Numerical problems. The main points consist of generating accurately com- 
plicated Bang-Singular sequences solutions of the Maximum Principle and to 
prove the convergence of the continuation problems. In this setting, the prob- 
lem is to initialize the shooting equation using the continuation method. 

Improving experimental results. The preliminary experimental figure IT51 
shows the problem of the magnetic fields inhomogeneities, which have therefore 
to be taken into account in the model. In this case, the geometric techniques 
can be used as a first step to initialize a purely iterative numerical approach 
such as the GRAPE algorithm [HI HH1 [2H 123 . In this setting, the geometric 
solution provides an efficient initial solution and gives the physical limit of the 
contrast problem that can be reached. The GRAPE algorithm is then able to 
solve the simultaneous optimal control of a large number of spins of an inho- 
mogeneous ensemble. In the example where the goal is to saturate the spins 
in deoxygenated blood, while maximizing the final magnetization of oxygenated 
blood, the numerically optimized pulse achieved about 93% and 70% of the con- 
trast found by the optimal geometric solution for the ideal and the real cases. 
Note that this problem with magnetic field inhomogeneities is related to the 
controllability analysis of [5]. 
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